title = "Venezuela")
ggarrange(F_Turkey, F_Venezuela, nrow = 2) + +
theme(axis.text=element_text(size=14),
axis.title=element_text(size=16,face="bold"))
ggarrange(F_Turkey, F_Venezuela, nrow = 2) +
theme(axis.text=element_text(size=14),
axis.title=element_text(size=16,face="bold"))
ggsave("outputs/Democratization_operationalization.pdf",  width = 35,
height = 30,
units = "cm",
dpi = 1200)
#### Subset Data ####
data_figure <- vdem_data %>%
filter(country_name %in% c("Turkey", "Venezuela"))
data_figure_turkey <- data_figure %>%
filter(year >=1900) %>%
relocate(country_name, year, dem_ep, dem, dem_bin_EDI) %>%
filter(country_name == "Turkey")
F_Turkey <- ggplot(data_figure_turkey, aes(x = year)) +
geom_line(aes(y = v2x_polyarchy)) +
geom_line(data =data_figure_turkey[data_figure_turkey$year >= 1907  & data_figure_turkey$year <=1909,],
aes(x = year, y = v2x_polyarchy), color = "red", size =1) +
geom_line(data = data_figure_turkey[data_figure_turkey$year >= 1945  & data_figure_turkey$year <=1951,],
aes(x = year, y = v2x_polyarchy), color = "red", size =1) +
geom_line(data = data_figure_turkey[data_figure_turkey$year >= 1961  & data_figure_turkey$year <=1967,],
aes(x = year, y = v2x_polyarchy), color = "red", size =1) +
geom_line(data = data_figure_turkey[data_figure_turkey$year >= 1983  & data_figure_turkey$year <=1988,],
aes(x = year, y = v2x_polyarchy), color = "red", size =1) +
geom_hline(yintercept = 0.5, linetype = "dashed") +
geom_label(
label="liberalizing\n autocracy",
x=1907,
y=0.2,
label.padding = unit(0.55, "lines"), # Rectangle size around label
label.size = 0.35,
color = "black",
fill="#69b3a2"
) +
geom_label(
label="liberalizing\n autocracy",
x=1948,
y=0.5,
label.padding = unit(0.55, "lines"), # Rectangle size around label
label.size = 0.35,
color = "black",
fill="#69b3a2"
) +
geom_label(
label="democratic\n transition",
x=1964,
y=0.65,
label.padding = unit(0.55, "lines"), # Rectangle size around label
label.size = 0.35,
color = "black",
fill="#AA4371"
) +
geom_label(
label="democratic\n transition",
x=1990,
y=0.25,
label.padding = unit(0.55, "lines"), # Rectangle size around label
label.size = 0.35,
color = "black",
fill="#AA4371"
) +
geom_point(aes(y=dem), color = "#453781FF") +
annotate(
"text", label = "Binary democracy\n measure (Step 1)",
x=1960,
y=0.90, size = 5, colour = "#453781FF" ) +
geom_point(aes(y=dem_bin_EDI+0.03), color = "#238A8DFF", alpha =0.5, shape = 2) +
annotate(
"text", label = "Binary democracy\n measure (Step 2)",
x=1925,
y=0.25, size = 5, colour = "#238A8DFF" ) +
theme_pubr() +
ylim(0,1.05) +
labs(y = "Electoral Democracy Index",
x = "Year",
color = "",
title = "Turkey")  +
theme(axis.text=element_text(size=14),
axis.title=element_text(size=16,face="bold"))
## Venezuela ##
data_figure_venezuela <- data_figure %>%
filter(year >=1900) %>%
relocate(country_name, year, dem_ep, dem, dem_bin_EDI) %>%
filter(country_name == "Venezuela")
F_Venezuela  <- ggplot(data_figure_venezuela, aes(x = year)) +
geom_line(aes(y = v2x_polyarchy)) +
geom_line(data =data_figure_venezuela[data_figure_venezuela$year >= 1945  & data_figure_venezuela$year <=1948,],
aes(x = year, y = v2x_polyarchy), color = "red", size =1) +
geom_line(data = data_figure_venezuela[data_figure_venezuela$year >= 1957  & data_figure_venezuela$year <=1960,],
aes(x = year, y = v2x_polyarchy), color = "red", size =1) +
geom_hline(yintercept = 0.5, linetype = "dashed") +
geom_label(
label="liberalizing\n autocracy",
x=1948,
y=0.5,
label.padding = unit(0.55, "lines"), # Rectangle size around label
label.size = 0.35,
color = "black",
fill="#69b3a2"
) +
geom_label(
label="democratic\n transition",
x=1968,
y=0.5,
label.padding = unit(0.55, "lines"), # Rectangle size around label
label.size = 0.35,
color = "black",
fill="#AA4371"
) +
geom_point(aes(y=dem), color = "#453781FF") +
annotate(
"text", label = "Binary democracy\n measure (Step 1)",
x=1960,
y=0.90, size = 5, colour = "#453781FF" ) +
geom_point(aes(y=dem_bin_EDI+0.03), color = "#238A8DFF", alpha =0.5, shape = 2) +
annotate(
"text", label = "Binary democracy\n measure (Step 2)",
x=1925,
y=0.25, size = 5, colour = "#238A8DFF" ) +
theme_pubr() +
ylim(0,1.05) +
labs(y = "Electoral Democracy Index",
x = "Year",
color = "",
title = "Venezuela") +
theme(axis.text=element_text(size=14),
axis.title=element_text(size=16,face="bold"))
ggarrange(F_Turkey, F_Venezuela, nrow = 2)
ggsave("outputs/Democratization_operationalization.pdf",  width = 35,
height = 30,
units = "cm",
dpi = 1200)
#### "Reanalyzing the Link between Democracy and Economic Development" ####
# authors: "Pelke, Lars"
# date: 2023-05-22
# written under "R version 4.1.2 (2021-11-01)"
#### Preliminaries ####
library(plm)
library(MASS)
library(multiwayvcov)
library(ggeffects)
library(PanelMatch)
library(tidyverse)
library(ggpubr)
#run code 01_data_wrangling and before using this code
# set your working directory here #
setwd("C:/Users/ba72loko/Promotion/Paper/Autocratization and Economic Well-Being/calculations_Nov_2021")
# clear workspace
rm(list=ls())
#### Import Data ####
# Load VDem data
vdem <- readRDS("data/vdem_data.rds")
vdem <- vdem %>%
filter(year >=1896)
#### Difference in Difference Design #### compare Imai et al. 2020, package: https://github.com/insongkim/PanelMatch
vdem$year <- as.integer(vdem$year)
vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations
#### Define outcome variable ####
vdem$y <- vdem$gdppc_ppp_bcbt_mean_log*100
summary(vdem$y)
##### Democratization and Economic Development ######
#### Panel Match Models ####
vdem <- as.data.frame(vdem)
class(vdem)
## Before Refinement ##
PM.results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "country_id",
treatment = "dem", refinement.method = "none",
data = vdem, match.missing = FALSE,
covs.formula = ~ I(lag(y, 1:4)),
size.match = 10, qoi = "att" ,outcome.var = "y",
lead = 0:20, forbid.treatment.reversal = FALSE,
use.diagonal.variance.matrix = TRUE)
## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 4, time.id = "year", unit.id = "country_id",
treatment = "dem", refinement.method = "mahalanobis",
data = vdem, match.missing = TRUE,
covs.formula = ~ I(lag(v2x_polyarchy, 1:4)) + I(lag(Maddison2020_pop_mean_log, 1:4)) +
I(lag(v2svstterr, 1:4)) + I(lag(v2clrspct, 1:4)) +
I(lag(gdppc_ppp_bcbt_growth_rate_5avg, 1:4)) +
I(lag(y, 1:4)),
size.match = 10, qoi = "att" ,outcome.var = "y",
lead = 0:20, forbid.treatment.reversal = FALSE,
use.diagonal.variance.matrix = TRUE)
## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 4, time.id = "year", unit.id = "country_id",
treatment = "dem", refinement.method = "ps.match",
data = vdem, match.missing = TRUE,
covs.formula = ~ I(lag(v2x_polyarchy, 1:4)) + I(lag(Maddison2020_pop_mean_log, 1:4)) +
I(lag(v2svstterr, 1:4)) + I(lag(v2clrspct, 1:4)) +
I(lag(gdppc_ppp_bcbt_growth_rate_5avg, 1:4)) +
I(lag(y, 1:4)),
size.match = 10, qoi = "att" ,outcome.var = "y",
lead = 0:20, forbid.treatment.reversal = FALSE,
use.diagonal.variance.matrix = TRUE)
## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 4, time.id = "year", unit.id = "country_id",
treatment = "dem", refinement.method = "ps.weight",
data = vdem, match.missing = TRUE,
covs.formula = ~ I(lag(v2x_polyarchy, 1:4)) + I(lag(Maddison2020_pop_mean_log, 1:4)) +
I(lag(v2svstterr, 1:4)) + I(lag(v2clrspct, 1:4)) +
I(lag(gdppc_ppp_bcbt_growth_rate_5avg, 1:4)) +
I(lag(y, 1:4)),
size.match = 10, qoi = "att" ,outcome.var = "y",
lead = 0:20, forbid.treatment.reversal = FALSE,
use.diagonal.variance.matrix = TRUE)
## CBPS Weighting ##
PM.results_cbpsweight <- PanelMatch(lag = 4, time.id = "year", unit.id = "country_id",
treatment = "dem", refinement.method = "CBPS.weight",
data = vdem, match.missing = TRUE,
covs.formula = ~ I(lag(v2x_polyarchy, 1:4)) + I(lag(Maddison2020_pop_mean_log, 1:4)) +
I(lag(v2svstterr, 1:4)) + I(lag(v2clrspct, 1:4)) +
I(lag(gdppc_ppp_bcbt_growth_rate_5avg, 1:4)) +
I(lag(y, 1:4)),
size.match = 10, qoi = "att" ,outcome.var = "y",
lead = 0:20, forbid.treatment.reversal = FALSE,
use.diagonal.variance.matrix = TRUE)
## CBPS Matching ##
PM.results_cbpsmatch <- PanelMatch(lag = 4, time.id = "year", unit.id = "country_id",
treatment = "dem", refinement.method = "CBPS.match",
data = vdem, match.missing = TRUE,
covs.formula = ~ I(lag(v2x_polyarchy, 1:4)) + I(lag(Maddison2020_pop_mean_log, 1:4)) +
I(lag(v2svstterr, 1:4)) + I(lag(v2clrspct, 1:4)) +
I(lag(gdppc_ppp_bcbt_growth_rate_5avg, 1:4)) +
I(lag(y, 1:4)),
size.match = 10, qoi = "att" ,outcome.var = "y",
lead = 0:20, forbid.treatment.reversal = FALSE,
use.diagonal.variance.matrix = TRUE)
######################################################
## Display different matched sets ##
mset <- PM.results_mahalanobis$att[9]
DisplayTreatment(unit.id = "country_id",
time.id = "year", legend.position = "none",
xlab = "year", ylab = "Country Code",
treatment = "dem", data = vdem,
matched.set = mset, # this way we highlight the particular set
show.set.only = TRUE)
summary(PM.results_mahalanobis$att)
summary(PM.results_psmatch$att)
summary(PM.results_psweight$att)
summary(PM.results_cbpsmatch$att)
summary(PM.results_cbpsweight$att)
par(mfrow=c(1,5))
plot(PM.results_mahalanobis$att, main ="Distribution of Matched Set Sizes\n (Mahalanobis Matching)")
plot(PM.results_psmatch$att,main = "Distribution of Matched Set Sizes\n (Propenstiy ScoreMatching)")
plot(PM.results_psweight$att,main = "Distribution of Matched Set Sizes\n (Propenstiy Score Weighting)")
plot(PM.results_cbpsmatch$att,main = "Distribution of Matched Set Sizes\n (Covariate Propensity Score Matching)")
plot(PM.results_cbpsweight$att,main = "Distribution of Matched Set Sizes\n (Covariate Propenstiy Score Weighting)")
dev.copy(pdf,'outputs/Step4/Figure_Frequency_Distribution_log.pdf', width = 10, height = 5)
dev.off()
#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####
## Get covariate balance ##
df_none <- get_covariate_balance(PM.results.none$att,
data = vdem,
covariates = c("v2x_polyarchy", "Maddison2020_pop_mean_log",  "v2svstterr", "v2clrspct",  "gdppc_ppp_bcbt_growth_rate_5avg", "y"),
plot = FALSE)
df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
data = vdem,
covariates = c("v2x_polyarchy", "Maddison2020_pop_mean_log",  "v2svstterr", "v2clrspct", "gdppc_ppp_bcbt_growth_rate_5avg", "y"),
plot = FALSE)
df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
data = vdem,
covariates = c("v2x_polyarchy", "Maddison2020_pop_mean_log",  "v2svstterr", "v2clrspct",  "gdppc_ppp_bcbt_growth_rate_5avg", "y"),
plot = FALSE)
df_psweight <- get_covariate_balance(PM.results_psweight$att,
data = vdem,
covariates = c("v2x_polyarchy", "Maddison2020_pop_mean_log",  "v2svstterr", "v2clrspct", "gdppc_ppp_bcbt_growth_rate_5avg", "y"),
plot = FALSE)
df_cbpsmatch <- get_covariate_balance(PM.results_cbpsmatch$att,
data = vdem,
covariates = c("v2x_polyarchy", "Maddison2020_pop_mean_log",  "v2svstterr", "v2clrspct",  "gdppc_ppp_bcbt_growth_rate_5avg", "y"),
plot = FALSE)
df_cbpsweight <- get_covariate_balance(PM.results_cbpsweight$att,
data = vdem,
covariates = c("v2x_polyarchy", "Maddison2020_pop_mean_log",  "v2svstterr", "v2clrspct", "gdppc_ppp_bcbt_growth_rate_5avg", "y"),
plot = FALSE)
#Tidying outputs #
df_none <- as.data.frame(df_none)
df_none <- df_none %>%
rownames_to_column("year") %>%
pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
mutate(year = case_when(year =="t_4" ~ -4,
year =="t_3" ~ -3,
year =="t_2" ~ -2,
year =="t_1" ~ -1,
year == "t_0" ~ 0))
F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
geom_line() +
geom_hline(yintercept=0, color = "black") +
labs(x = "", y = "",
title = "Before Refinement") +
ylim(-2, 5) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
theme_pubr()
df_mahalanobis <- as.data.frame(df_mahalanobis)
df_mahalanobis <- df_mahalanobis %>%
rownames_to_column("year") %>%
pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
mutate(year = case_when(year =="t_4" ~ -4,
year =="t_3" ~ -3,
year =="t_2" ~ -2,
year =="t_1" ~ -1,
year == "t_0" ~ 0))
F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
geom_line() +
geom_hline(yintercept=0, color = "black") +
labs(x = "", y = "",
title = "Mahanlanobis Distance\n Matching") +
ylim(-2, 5) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
theme_pubr()
df_psmatch <- as.data.frame(df_psmatch)
df_psmatch <- df_psmatch %>%
rownames_to_column("year") %>%
pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
mutate(year = case_when(year =="t_4" ~ -4,
year =="t_3" ~ -3,
year =="t_2" ~ -2,
year =="t_1" ~ -1,
year == "t_0" ~ 0))
F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
geom_line() +
geom_hline(yintercept=0, color = "black") +
labs(x = "", y = "",
title = "Propensity Score\n Matching") +
ylim(-2, 5) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
theme_pubr()
df_psweight <- as.data.frame(df_psweight)
df_psweight <- df_psweight %>%
rownames_to_column("year") %>%
pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
mutate(year = case_when(year =="t_4" ~ -4,
year =="t_3" ~ -3,
year =="t_2" ~ -2,
year =="t_1" ~ -1,
year == "t_0" ~ 0))
F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
geom_line() +
geom_hline(yintercept=0, color = "black") +
labs(x = "", y = "",
title = "Propensity Score\n Weighting") +
ylim(-2, 5) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
theme_pubr()
df_cbpsmatch <- as.data.frame(df_cbpsmatch)
df_cbpsmatch <- df_cbpsmatch %>%
rownames_to_column("year") %>%
pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
mutate(year = case_when(year =="t_4" ~ -4,
year =="t_3" ~ -3,
year =="t_2" ~ -2,
year =="t_1" ~ -1,
year == "t_0" ~ 0))
F5 <- ggplot(df_cbpsmatch, aes(x = year, y = sd_covariate, color = covariate)) +
geom_line() +
geom_hline(yintercept=0, color = "black") +
labs(x = "", y = "",
title = "Covariate Balancing Propensity Score\n Matching") +
ylim(-2, 5) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
theme_pubr()
df_cbpsweight <- as.data.frame(df_cbpsweight)
df_cbpsweight <- df_cbpsweight %>%
rownames_to_column("year") %>%
pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
mutate(year = case_when(year =="t_4" ~ -4,
year =="t_3" ~ -3,
year =="t_2" ~ -2,
year =="t_1" ~ -1,
year == "t_0" ~ 0))
F6 <- ggplot(df_cbpsweight, aes(x = year, y = sd_covariate, color = covariate)) +
geom_line() +
geom_hline(yintercept=0, color = "black") +
labs(x = "", y = "",
title = "Covariate Balanacing Propensity Score\n Weighting") +
ylim(-2, 5) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
theme_pubr()
Figure1 <- ggarrange(F1, F2, F3, F4, F5, F6,
ncol = 3, nrow = 2,
common.legend = TRUE)
annotate_figure(Figure1,
left = text_grob("Standardized Mean Differences\n for Democratization", rot = 90))
ggsave("outputs/Step4/Figure_Covariate_Balance_Autocratization.pdf", height = 15, width = 35, units= c("cm"))
#### Empirical Findings of Democratization Treatment ####
## Model 1 ##
M1 <- PanelEstimate(sets = PM.results_mahalanobis, data = vdem)
summary(M1)
par(mfrow=c(1,1))
plot(M1)
## Model 2 ##
M2 <- PanelEstimate(sets = PM.results_psmatch, data = vdem)
summary(M2)
## Model 3 ##
M3 <- PanelEstimate(sets = PM.results_psweight, data = vdem)
summary(M3)
## Model 4 ##
M4 <- PanelEstimate(sets = PM.results_cbpsweight, data = vdem)
summary(M4)
## Model 5 ##
M5 <- PanelEstimate(sets = PM.results_cbpsmatch, data = vdem)
summary(M5)
#### Plot the results
par(mfrow=c(2,3))
plot(M1, main = "Mahalanobis Matchting", ylab= "Estimated Effect of Democratization", xlab = "", ylim = c(-5, 50))
plot(M2, main = "Propensity Score Matching", ylab= "", xlab = "", ylim =  c(-5, 50))
plot(M3, main = "Propensity Score Weighting", ylab= "", xlab = "", ylim =  c(-5, 50))
plot(M4, main = "CBPS Weighting", ylab= "Estimated Effect of Democratization", xlab = "", ylim =  c(-5, 50))
plot(M5, main = "CBPS Matching", ylab= "", xlab = "", ylim =  c(-5, 50))
dev.copy(pdf,'outputs/Step4/A_Figures_Democratization_Estimates_log.pdf', width = 10, height =6)
dev.off()
## Plotting the results in Paper Design ##
M1_df <- data.frame(M1$estimates, M1$standard.error, M1$lead)
M1_df$min95 <-   M1$estimates- 1.96*M1$standard.error
M1_df$max95 <-   M1$estimates+ 1.96*M1$standard.error
View(M1_df)
#### "Reanalyzing the Link between Democracy and Economic Development" ####
# authors: "Pelke, Lars"
# date: 2021-11-26
# written under "R version 4.1.2 (2021-11-01)"
#### Preliminaries ####
R.version$version.string
# clear workspace
rm(list=ls())
# set working directory
setwd("C:/Users/ba72loko/Promotion/Paper/Autocratization and Economic Well-Being/calculations_final")
# loading packages
library(countrycode)
library(tidyverse)
library(viridis)
library(ggpubr)
library(texreg)
library(sjPlot)
library(ggeffects)
library(readstata13)
#library(vdemdata)
library(imputeTS)
library(stargazer)
library(haven)
#### Import Data ####
# Load VDem data
vdem <- readRDS("data/vdem_10/Country_Year_V-Dem_Full+others_R_v10/V-Dem-CY-Full+Others-v10.rds")
#### Regime duration ####
vdem <- vdem %>%
mutate(autocracy= ifelse(v2x_regime<2, 1, 0),
autocracy= ifelse(is.na(v2x_regime), NA, autocracy))
cumsum.na <- function(x) {
x[which(is.na(x))] <- 0
return(cumsum(x))
}
vdem_auto <- vdem %>%
group_by(country_id) %>%
mutate(democracy = ifelse(autocracy==1, 0, 1),
democracy = ifelse(is.na(v2x_regime), NA, democracy)) %>% # get binary democracy measure for constructing years under democracy
mutate(autocracy_years = ifelse(autocracy==1, cumsum.na(autocracy), NA)) %>%  # years under autocracy
mutate(democracy_years = ifelse(democracy==1, cumsum.na(democracy), NA)) %>% # years under democracy
ungroup() %>%
dplyr::select(country_id, year, autocracy_years, democracy_years)
vdem_auto$autocracy_years_na <- replace(vdem_auto$autocracy_years, is.na(vdem_auto$autocracy_years), 0)
vdem_auto$democracy_years_na <- replace(vdem_auto$democracy_years, is.na(vdem_auto$democracy_years), 0)
vdem_auto <- vdem_auto %>%
mutate(regime_duration = democracy_years_na,
regime_duration = regime_duration + autocracy_years_na)
vdem_auto <-vdem_auto %>%
dplyr::select(country_id, year, autocracy_years, democracy_years, regime_duration) %>%
mutate(regime_duration = ifelse(year<1900, NA, regime_duration))
vdem <- vdem %>%
left_join(vdem_auto, by = c("country_id", "year"))
vdem <- vdem %>%
distinct(country_name, year, .keep_all = TRUE)
vdem$cowcode[vdem$country_name == "Hong Kong"] <- 715
vdem$cowcode[vdem$country_name == "Palestine/West Bank"] <- 1020
#### Episodes of Regimes Dataset ####
devtools::install_github("vdeminstitute/ERT")
ert_data <- ERT::get_eps(
data = vdemdata::vdem,
start_incl = 0.05,
cum_incl = 0.1,
year_turn = 0.03,
cum_turn = 0.1,
tolerance = 5
)
#### "Reanalyzing the Link between Democracy and Economic Development" ####
# authors: "Pelke, Lars"
# date: 2023-05-22
# written under "R version 4.1.2 (2021-11-01)"
#### Preliminaries ####
library(plm)
library(MASS)
library(multiwayvcov)
library(ggeffects)
library(PanelMatch)
library(tidyverse)
library(ggpubr)
#run code 01_data_wrangling and before using this code
# set your working directory here #
setwd("C:/Users/ba72loko/Promotion/Paper/Autocratization and Economic Well-Being/calculations_final")
# clear workspace
rm(list=ls())
#### Import Data ####
# Load VDem data
vdem <- readRDS("data/vdem_data.rds")
